#include "Tools.h"

//calculates the norm of a vector
//=============================================
//IN:
//pointer with components x, y and z.
//OUT:
//the norm
double norm(const double *n);

double setScalarField(double *xyz,int option);

double norm(const double *n){
	return sqrt( pow(n[0],2) + pow(n[1],2) + pow(n[2],2) );
}

double F_area(const double *ptn1, const double *ptn2, const double *ptn3){

	double xyz[3][3], n[3];
	for (int i=0;i<3;i++){
		xyz[0][i] = ptn1[i];
		xyz[1][i] = ptn2[i];
		xyz[2][i] = ptn3[i];
	}
	double a[3] = { xyz[1][0]-xyz[0][0], xyz[1][1]-xyz[0][1], xyz[1][2]-xyz[0][2] } ;
	double b[3] = { xyz[2][0]-xyz[0][0], xyz[2][1]-xyz[0][1], xyz[2][2]-xyz[0][2] } ;
	computeCrossProduct(a,b,n);

	return .5*norm(n);
}

double R_Volume(const double *ptn1, const double *ptn2, const double *ptn3, const double *ptn4){

	double x[4] = {ptn1[0],ptn2[0],ptn3[0],ptn4[0]};
	double y[4] = {ptn1[1],ptn2[1],ptn3[1],ptn4[1]};
	double z[4] = {ptn1[2],ptn2[2],ptn3[2],ptn4[2]};
	double a[3] = {x[1]-x[0],y[1]-y[0],z[1]-z[0]} ;
	double b[3] = {x[2]-x[0],y[2]-y[0],z[2]-z[0]} ;
	double c[3] = {x[3]-x[0],y[3]-y[0],z[3]-z[0]} ;
	double normal[3] ;
	computeCrossProduct(a,b,normal);
	double value = .1666666666667*computeDotProduct(normal,c);
	return fabs(value);
}

void generateScalarField(pMesh theMesh, DataTransfer data_transfer, int option){

	//cout<<"\n\n#Printing "<<": \n"<<endl;
	long double pi = 3.141592653589793;
	//loop over vertices
	VIter vit = M_vertexIter(theMesh);
	while (pEntity vertex = VIter_next(vit)){

		int id = EN_id(vertex);
		//cout<<"Node "<<id<<": ";

		//get coordinates
		double xyz[3] = {0,0,0};
		V_coord(vertex, xyz);

		double value = setScalarField(xyz,option);
		//cout<<"-> FIELD: "<<value<<endl;
		//store scalar acoording to the function above
		data_transfer.setNodeValue(id,value);
		
		//cout<<data_transfer.getNodeValue(id)<<endl;

	}
	VIter_delete(vit);


}

void calculate_Error(pMesh theMesh, DataTransfer data_transfer, int option) {
  
  //variables 
  int IDmax;
  double maxError=0; double summ = 0; double summ2 = 0;
  
  VIter vit = M_vertexIter(theMesh);
  while (pEntity vertex = VIter_next(vit)){
        
    int id = EN_id(vertex);
    //get coordinates
    double xyz[3] = {0,0,0};
    V_coord(vertex, xyz);
    
    double realValue = setScalarField(xyz,option);
    double interpValue = data_transfer.getInterpolatedValues(id);
    
    //node error
    double error = (realValue - interpValue)/realValue;
    error= fabs(error);
    
    //maximum norm
    if (error > maxError) maxError = error;
    
    double error2 = realValue - interpValue;
    summ = summ + pow(error2,2);
    summ2= summ2 + pow(realValue,2);    
  }
  
  //calculating the 2-norm
  summ=sqrt(summ); summ2= sqrt(summ2);
  
  //printing
  cout<<"Max. Node Error "<<maxError<<"%"<<endl;
  cout<<"Global error (2-norm): "<<summ/summ2<<"%"<<endl;  
  
  VIter_delete(vit);
}

double setScalarField(double *xyz,int option){
  
  //create function
  double value;
  switch (option){
    case 1: value = 10;break;
    case 2: value = 5*xyz[0] + 10; break;
    
    /*case 3: {
      
      if(xyz[1] < 5) value = 20;
      else value = 10;
      break;
    }*/
    case 3: value = pow(xyz[1],2) + 20; break;
    case 4: value = pow(xyz[0],3) + xyz[0] + 4;break;
  }
  
  return value;
}
  